#pragma once
#include "MatrixV.h"
#include "function.h"
#include <lapacke.h>
using namespace std;
/// (I-theta*k*A)U^{n+1} = (I+(1-theta)*k*A)U^n
class theta_method
{
private:
    /// 几个参数的关系为h = 1/(m+1), r = (k*nu)/(h*h)
    int m;
    double theta, r, h, nu, k;
    vector<vector<double>> M1, M2;
public:
    theta_method(double _theta, double _r, double _h):theta(_theta), r(_r), h(_h)
	{
	    m = (int)1/h - 1;
	    nu = 1;
	    k = h*h*r/nu;
	    MatrixV M(r, nu, h);
	    vector<vector<double>> A = M.A(), I = M.I();
	    M1 = I - theta*k*A;
	    M2 = I + (1 - theta)*k*A;
	};
    /// 根据初始条件函数创建 U^0
    vector<double> U0()
	{
	    vector<double> U0;
	    for (int i = 0; i < m; i++)
	    {
		U0.push_back(u0((i+1)*h));
	    }
	    return U0;
	};
    /// 求解器,输入时间步数,保存求解结果到 U
    vector<double> Solver(int t)
	{
	    vector<double> UN(m), MU;
	    vector<vector<double>> U;
	    U.push_back(U0());
	    for (int i = 0; i < t; i++)
	    {
		MU = M2*U[i];/// 问题变为 M1*UN = MU -> AU = F;
		UN = solveU(M1,MU);
                // for (int j = 0; j < m; j++)
		// {
		//     F[j] = MU[j];
		//     for (int l = 0; l < m; l++)
		//     {
		// 	A[j+l*m] = M1[j][l];
		//     }
		// }
		// lapack_int info,m0,n0,lda,ldb,nrhs;
		// m0 = m;
		// n0 = m;
		// nrhs = 1;
		// lda = m;
		// ldb = m;
		// info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m0,n0,nrhs,A,lda,F,ldb);
		// for(int i = 0; i < m; i++)
		// {
		//     UN[i] = F[i];
		// }
		// print(UN);
		U.push_back(UN);
	    }
	    return U[t];
	};
};
